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One may predict a quasicrystal structure starting from 
electrons and quantum mechanics, as approximated by inter- 
atomic pair potentials calibrated with ab-initio total-energy 
calculations, combined with the experimentally known com- 
position and lattice constants. Here we report our progress on 
the "basic Ni" decagonal phase d(Al7oNi2iCog. Atomic con- 
figurations are represented as decorations of (possibly) ran- 
dom tilings. Our method was Monte Carlo simulation using 
both lattice-gas hops by atoms and tile-flip rearrangements, 
eventually followed by molecular dynamics and relaxation of 
the atom positions. Initially allowing the greatest freedom of 
atom positions, we observed nearly deterministic structural 
rules and enforced these as constraints involving larger tiles; 
this procedure was repeated at the next level of modeling. 
In crude and preliminary form, the effective Hamiltonian for 
tile-tile interactions is known, which is needed for further sim- 
ulations to infer the long-range order. Our atomic arrange- 
ments in the 20 A decagonal cluster are compared with three 
structure models based on recent experiments. Q 



I. INTRODUCTION 

This paper reports on early results of a program to 
determine the structure of a quasicrystal - specifically 
d(AlNiCo) - strictly from energetics, without use of 
diffraction data. Almost all of the work was performed 
by M.M. and M.W., with their coauthors, and is reported 
in Ref. [Q, but the present account is tilted towards the 
particular interests of C.L.H. 

We have focused on the thermodynamically sta- 
ble and highly ordered "Basic Nickel" subphase of 
decagonal d(AlNiCo) of composition approximately 
Alo.70Nio.21Coo.09- Our ultimate objective is to predict 
the structure on the basis of total energy; though the 
gross aspects of the structure are clear from diffraction, 
there are many ambiguous or controversial details, and 
there is no understanding for the quasicrystal's special 
stability and simplicity at this composition. In particu- 
lar, we wish to identify the sites of the transition metal 
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("TM") components Ni and Co, which are indistinguishi- 
ble to (ordinary) X-rays or electrons j|-|§ ■ 

In a sense, the most difficult problem of crystal chem- 
istry is to predict which structure a given composition 
will adopt, even if one has an exact and tractable Hamil- 
tonian for the total energy (e.g. by pair potentials) - 
for this necessarily involves a comparison with an in- 
finity of possible structures. The only irreproachable 
method is a mathematical proof, which is feasible in rare 
cases with short-range interactions, e.g. close-packing 
of hard spheres, or the two-dimensional "binary-tiling" 
quasicrystal (a toy model) A naive, brute- force ap- 
proach would be to cool (in simulation) from the melt, 
and see which structures emerge; but this is prone to 
fail for a complicated material, since the accessible times 
are so short that it will get stuck in a glassy disordered 
configuration. Even when a quasicrystal emerges, as 
with Dzugutov's toy potentials [[10| , the system may have 
found a merely metastable ordered state for kinetic rea- 
sons: indeed, the stable phase in this case turned out to 
be a simple bcc packing fl~ij|. 



II. DECORATION MODELS 

We represent the quasicrystal structure as a decoration 
of disjoint tiles. There are several reasons for doing so, 
rather than as a irrational cut through a five-dimensional 
hypercrystal. Based in real space, this representation is 
easier to visualize and somewhat more tractable, tech- 
nically, than the hypercrystal approach. For example, 
shifts of the atoms from ideal (tiling-vertex) sites are 
parametrized by just a finite set of real numbers, like the 
coordinates in a crystal unit cell. A decoration descrip- 
tion - particularly in the decagonal case - naturally lends 
itself to a hierarchy of supertilings, which are presumably 
involved in the modulations which distinguish some of 
the subphases in the d(AlNiCo) phase diagram ||. Fur- 
thermore, most hyperspace structures which have been 
defined using discrete acceptance domains, e.g. Ref. [0, 
can be translated into tiling-decoration language. 

Decoration models are uniquely suited to a correct 
modeling of chemical disorder (including vacancies) in 
a quasicrystal. Of course, a hyperspace fit to Bragg data 
can incorporate mixed occupation of a site type, but it 
is highly implausible (from the viewpoint of the struc- 
tural energy) that the occupations of neighboring sites 
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vary independently of each other. If the occupations are 
strongly correlated, then the mean structure may differ 
significantly from any particular real structure; this mat- 
ters e.g. because the electronic structure is quite sensitive 
to the TM-TM contacts in d(AlNiCo) @. The tiling 
decoration is a very convenient way to account for such 
occupation correlations, by ascribing the randomness to 
the tiles and not the decoration. The ultimate criterion 
of a decorated model - fixing e.g. the appropriate size of 
tile - is that the ensemble of allowed tilings corresponds 
one-to-one with the ensemble of low-energy arrangements 
in an atomistic model. 

A decoration description replaces the independent real 
coordinates of many atoms by tiling degrees of freedom, 
which are discrete and many times fewer. When the dec- 
oration is deterministic and disjoint (each atom bound to 
a particular tile or tile-related geometric object), then the 
atoms' interaction energies can be expressed as a function 
solely of the tile configuration, called the "tile Hamilto- 
nian" Htn c - This facilitates faster Monte Carlo simula- 
tions, on the comparatively few tile degrees of freedom, 
in order (i) to reveal supertilings in the process of dis- 
covering the atomic structure, as in this work; (ii) to 
predict diffuse scattering, or (iii) eventually to measure 
the phason elastic constants (which govern the diffuse 
wings surrounding Bragg peaks.) Armed with a reliable 
7itiic, one can also simulate three-dimensional samples 
to obtain quantitative estimates of the perfection of the 
long-range order and the strength of diffuse diffraction. 

This approach may also be the easiest way to resolve, 
for a particular stable quasicrystal, whether or not it 
is entropically stabilized. It is stabilized by energy fa- 
voring a quasiperiodic state if 7itiic happens to imple- 
ment the Penrose 'matching rules'. That almost occurs 
in d(AlCuCo), as modeled with potentials like those used 
in the present work [ p^[ . 

Our challenge, then, is to find the decoration rule and 
the tile-packing constraints without any bias, apart from 
the above-mentioned assumptions about the lattice con- 
stants. A caveat to keep in mind is that the results will 
probably be quite sensitive to the exact composition and 
number density. 

III. INPUTS: POTENTIALS AND CONSTRAINTS 

The main input to our calculations is a set of atom- 
atom pair potentials Vy(r) for species i and j, six dis- 
tinct functions in the case of a ternary alloy. They are 
given by Moriarty's "generalized pseudopotential theory" 
(GPT) p^-[l^], a systematic expansion of the total en- 
ergy as a sum Eq + Ex + E% + . . ., where the E n term 
depends on n atoms, but we use only the n = 2 terms. 
The first-neighbor well of the raw Vtmtm is unphysi- 
cally deep (which would be canceled by E3 in the sys- 
tematic theory). Therefore our Vtmtm potentials were 
empirically modified, by fitting a short-range repulsive 



correction so as to match all the forces in a full density- 
functional calculation on a small (50 atom) approximant 
of d( AINiCo) |f6| . There is strong support that our Vy (r) 
have the correct r dependence and relative strength, since 
they predict the Al-Co and Al-Ni phase diagrams pretty 
well |l5|] (as a function of concentration at T = 0); how- 
ever the ternary Al-Ni-Co phase diagram was not at- 
tempted. On the basis of the simulated melting tem- 
perature and phonon spectra, we do suspect that our 
potentials are ~30% stronger than the reality. 

Like the similar Al-TM potentials of Phillips, Zou, and 
Carlson, used earlier in Rcfs. jl3| and p7| , the GPT po- 
tentials depend implicitly on the net valence electron den- 
sity, and exhibit strong "Friedel oscillations" as a func- 
tion of distance (tails decaying as cos(2fcj?r + S)/r 3 ). For 
example, Vaitm(^) and Vtmtm (r) show four prominent 
minima in the range r < 8.5 A; such minima are quite 
important in deciding the structure [ p~7]JT8| . 

Notice that Vaiai(^) is practically zero after its hard 
core radius ~ 2.9 A : there is not even a nearest-neighbor 
well, just a shoulder on the tail of the hard-core. The 
Vaitm potentials have very deep nearest-neighbor wells: 
this creates the illusion of a TM-TM repulsion, but re- 
ally TM-TM neighbors are avoided mainly to make more 
room for Al in the TM coordination shells. However, 
at roughly 30% TM content, some TM-TM contacts are 
unavoidable; they are all Ni-Ni, since Vaico has a deeper 
well than Vain; ■ Finally, the Vtmtm potentials have their 
deepest well at second-neighbor (~ 4.2A) distances. The 
consequence of all this is that the TM atoms form a rather 
rigid and somewhat uniformly spaced network , while 
Al atoms move rather freely to follow the potential wells 
or troughs created by the TM arrangement. (See Sec. |v| 
for more discussion). 

Our use of these potentials imposes some limitations. 
We omitted the n = term, which depends only on the 
valence electron density, and is the largest contribution 
to a metal's total energy. Hence we can make valid com- 
parisons only between structures with (practically) the 
same valence electron density. We cannot meaningfully 
predict the lattice constants. Nor can we even expect 
the quasicrystal phase to be globally stable according to 
our potentials, since they might spuriously favor phase- 
separation. 

In view of the above-mentioned pitfalls for a brute- 
force approach, we adopted as a second input the ex- 
perimental lattice constants [Q: the stacking period is 
c = 4.08 A and the quasilattice constant in the decago- 
nal plane is oq — 2 Ah A. We only seek the lowest energy 
among arrangements on this framework. This highly con- 
straining assumption still permits a vast ensemble of pos- 
sible structures. 
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IV. METHODS AND RESULTS 

Our methods are a mix of tile-flip Monte Carlo, atom- 
hopping Monte Carlo, relaxation of atomic positions, 
and molecular dynamics. Our procedure is first to dis- 
cover the favorable low-energy motifs through Monte 
Carlo annealing, then to remove unnecessary degrees of 
freedom, and repeat, producing successively more con- 
strained models. Having fewer degrees of freedom the 
latter are much faster to simulate at low temperatures. 
At the end, we can investigate the effects of letting atoms 
depart from ideal sites. 

The initial stage of our exploration stacks two inde- 
pendent small (edge ao) Penrose rhombic tilings, in a 
vertical space of one lattice constant c = 4.08 A and 
the only allowed atom sites are on vertices. A manage- 
able size was 50 atoms on 72 candidate sites, initially on 
a good Pcnrosc-tiling approximant; periodic boundary 
conditions are used in all directions (and at all stages 
of exploration). The atoms - initially chosen to approx- 
imate the experimental density and composition of the 
Ni-rich "Basic Ni" phase - hop as a lattice gas on these 
sites. This allows sufficient freedom for the atoms, if they 
"want" , to adopt any of the decagonal structure models 
- with stacking period c = 4 A - that were hypothe- 
sized at one time or another 20 2W. The Monte Carlo 



moves permit swaps of the species between two nearby 
sites ("vacant" is a special case of species!), as well as 
"tile flips" which reshuffle the three rhombi in a fat or 
thin hexagon in the same layer. 

After this model is slowly cooled to zero temperature, if 
the initial composition was rightly chosen, one obtains a 
one-layer Hexagon-Boat-Star (HBS) tiling of edge length 
ao, with a two-layer decoration in which the allowed sites 
lie over vertices of the Penrose rhombi (into which HBS 
tiles may be decomposed) or an additional site in each 
Fat rhombus. This decoration places Al atoms over the 
HBS vertices (even and odd vertices alternate between 
even and odd layers). Isolated Co atoms sit over each 
Boat and Star tile center, ringed by possible Al sites; in 
the Star at most two Al can be present, out of five ideal 
sites in this ring. Particularly characteristic are the ver- 
tical zigzag chains, with Ni in each layer (appearing as 
NiNi pairs in Fig. yj), over the interior of every Hexagon 
tile. However, NiNi may be replaced by CoAl in some 
places. Decorating a tiling with the same H/B/S ratios 
as the quasiperiodic Penrose tiling yields an ideal compo- 
sition Alo.700Nio.207Coo.093 and atomic volume 14.16 A 3 . 
This composition coincides with the experimental Basic 
Ni phase, but the number density of atoms is at least 5% 
greater than in experiment. 

For the next stage of modeling, to discover larger-scale 
regularities, the small HBS tiles which emerged from the 
initial stage arc elevated to fundamental objects, either 
with a rigid decoration or with some atoms fixed and 
others forming a lattice gas. The allowed tile flips are 
reshuffiings of fat hexagons of the underlying rhombus 



tiling, provided the result is a valid HBS tiling and con- 
serves the atom content. Additionally the Al pair inside 
the Star tile can rotate among five allowed orientations. 
From this simulation, it emerges that the Hexagon tiles 
containing Ni chains only touch tip-to-tip, so that the 
angles relating them are multiples of 72°. In fact, (see 
Fig. |l|) the long axes of the Hexagons form edges of an 
HBS supertiling with an inflated edge length r 2 ao, where 

r = (i + VE)/2. 



V. THE 20 A DECAGON CLUSTER 

We now compare our results with well-known struc- 
ture models, organizing the discussion around the famous 
20 A diameter decagon cluster, which is prominent in 
electron-microscope images and in most diffraction re- 
finements as well. (The ideal decagon diameter is actu- 
ally 2t 3 Oo.) Such decagons indeed appear in our struc- 
tures [^]; Fig. [j] is centered on one of them. Z-contrast 
images, in which the intensity is a direct projection of 
atoms weighted by squared atomic number, reveal the 
TM positions 10,0; the TM (and most Al) positions pro- 
posed by Ref. are practically the same as in our version 
of the cluster. 

For a more detailed comparison, we focus on three re- 
cent structural studies All of these, and our sim- 
ulation, agree on the following details of the projection: 
(1) the outer decagon (edge length r 2 ao) has Al on each 
vertex and a pair of TM atoms on each edge, which we 
identify as Ni. (2) A middle decagon (edge too) has TM 
atoms on each vertex, which we identify as Co. (3) An 
inner decagon (edge ao) has Al on each vertex. (4) In 
the center, the 10-fold symmetry is broken and a sort of 
isosceles triangle is observed, with one TM (we say Co) 
at the unique corner and pairs of TM (we say NiNi) at 
the base corners. The cluster is evident in Ref. Q on the 
right side of their Fig. 7, as a combination of a Boat tile 
+ 2 Hexagon tiles, just as in our Fig. [j]. 

Detail (4) is somewhat controversial, since some 
decagon images have non-triangular centers. Indeed, the 
density maps from the refinement of Ref. Q show six 
strongly TM sites at the center of many (but not all) of 
the 20 A decagons. We ascribe this to stacking flips be- 
tween one layer and the next layer (see Sec. 
projection, as in the right-side decagon of Fig. 

The simulations we described up till now implicitly 
assumed a strict c = 4.08 A periodicity in the stacking 
direction, neglecting (like so much other modeling) the 
fact that decagonal quasicrystals are three-dimensional. 
To obtain the configuration in Fig. ^|, we increased the 
periodicity to 2c, i.e. two HBS layers, which initially 
were the same tiling. First a tile-flip was made in one 
of the two HBS layers, then the tilings were decorated, 
annealed by MD, and relaxed to an energy minimum. 

In a fifth, controversial detail of the decagons, our 
model initially disagreed with experiments in which ev- 
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ery edge of the middle decagon shows a pair of Al atoms. 
These sites form vertical zigzag chains, but they can- 
not all be occupied since they are separated by essen- 
tially the interlayer spacing c/2 = 2.08 A. In our fixed- 
decoration model, then, these sites are occupied in only 
one atom layer, with separation c, so only one atom is 
visible on each edge. We found that, after molecular dy- 
namics (MD) at 1000K followed by relaxation, some Al 
neighbors of the Co atoms shifted to the middle-decagon 
edges, which now present Al doublets in projection, in 
agreement with experiment. 

To further understand what happens with these atoms, 
the time-average of the Al positions during the MD sim- 
ulation is shown in Fig. [|. Some 40% of Al atoms are 
rather delocalized, and would need a highly anisotropic 
Debye- Waller factor in a crystallographic fit. H] The 
vertical projection (Fig. ^|(a) shows, consistent with our 
remarks on the potentials (Sec. |l|), that each Co atom in 
the middle decagon is surrounded by a potential trough 
in which Al atoms appear almost free to roll like ball- 
bearings p^] . 

A slice along a vertical plane further clarifies the Al be- 
havior: a second type of trough extends vertically, with 
a zigzag shape, and in fact connects with the circular 
troughs. In our simulations with a cell 2c in the verti- 
cal direction, we found three Al atoms appearing in each 
zigzag trough (per 2c, i.e. per 4 atom layers). Notice 
that the z displacement of two of these atoms makes the 
layer puckered, as is already known from diffraction 
and implies a local cell-doubling along the c axis, similar 
to the Ali3Fe4 or AI13C04 decagonal approximants. 

The symmetry-breaking of the cluster interior was re- 
cently predicted from energies by a full ab-initio calcula- 
tion |p6| . They find, from a quite different starting point 
than ours, the identical arrangement of the 5 Al + 5 TM 
atoms (per atomic bilayer) found at the cluster's center. 
However, they rejected the possibility of Al doublets on 
the middle decagon edges without trying z relaxations, 
nor did they address Co/Ni ordering. We find it striking 
that our mere pair potentials sufficed to obtain the ex- 
perimentally indicated form of symmetry breaking in the 
cluster center, while using no experimental input whatso- 
ever regarding the positioning of Al relative to TM atoms. 

VI. TILE HAMILTONIAN AND 
THREE-DIMENSIONAL STACKING 

To model a large, three-dimensional sample using a 
tile Hamiltonian, wc cannot demand perfect periodicity 
in the z direction. (One reason is that the "entropic 
stabilization" explanation of quasicrystal order requires 
an extensive entropy.) So, we represent each bilayer of 
atoms by a distinct HBS tiling, and constrain adjacent 
tilings in the stack to differ by "stacking flips" , which are 
exactly the same reshufflings that constitute our Monte 
Carlo tile-flips. 



We also found (from simulations) that - among (small) 
HBS tilings in which H tiles form a network touching 
tip-to-tip - the in-plane tile Hamiltonian is dominated 
by interactions between "H(NiNi)" tiles (Hexagon tiles 
decorated by NiNi). Thus we suggest 

Htile = C 72^72 + E* N S {. (1) 

where ./V72 is the number of H(NiNi) tile pairs related by 
a 72° rotation about their common tip, and 7V s f is the 
total number of stacking flips. 

We fitted c 72 = 0.218eV [0 from the simulation with 
ideal sites; this is positive since the 72° relation of 
H(NiNi) tiles creates infertile NiNi neighbors, and thus 
reduces the number of (favorable) AIM pairs. The struc- 
ture minimizes N72 by arranging that in some places, an 
H(NiNi) tile has one tip which contacts no other H(NiNi), 
and sees CoAl (at a 144° angle) in place of NiNi. That 
may be accomplished in two ways. If we re-allow lattice- 
gas hopping on the interior sites of the (small) HBS tiles, 
and properly adjust the stoichiometry [pj, we obtain the 
same old rigid decoration, except that a few H(CoAl) 
tiles appear, as illustrated in Fig. [t[ On the other hand, 
if we retain the rigid NiNi decoration of small H tiles, the 
CoAl belongs to the interior of a Boat tile, and the large 
HBS tiling acquires a fourth kind of tile ("Bowtie"), as 
in Fig. 4(a) of Ref. §. 

Existing experiments can in principle measure E s f. 
The time-dependent local "phason" flips [27]], observed 
at T w 1200K by TEM through very thin samples of 
d(AlCuCo), can only occur by nucleating a "stacking 
flip" at one surface, which subsequently random-walks 
to the other surface. Video images |27j| display an in- 
termediate state for a noticeable fraction - say 10% - 
of the time, so that with a sample thickness estimated 
at 100 A (25 bilayers), one would very roughly esti- 
mate E s { « 0.3 eV in d(AlCuCo). It would be exciting 
if such observations could be analyzed quantitatively. In 
our simulation, with rigid atoms decorating ideal sites, 
E s f « 1.4 eV. However, relaxation (as in Fig. [2]) nearly 
cancels this energy cost, indeed E s { ss —0.2 eV in very 
preliminary results, which should not be compared with 
the experiments on a different quasicrystal. (Physically, 
E s f < probably implies a period 2c ~ 8 Amodulation.) 

VII. DISCUSSION 

We find it remarkable that, with quite sketchy ex- 
perimental input, our simulations appear competitive 
with single-crystal diffraction as a way of discovering 
this rather complex structure. This occurs despite the 
shortcomings of the pair potentials - our omission of 3- 
body terms, and the likelihood that the Vaitm we use 
is stronger than the real one. Even if the connection to 
the microscopies is not quantitative, one may still obtain 
the large-scale order quite well: that depends mostly on 
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the tile Hamiltonian (Q) having the proper form, which is 
more robust than the numerical values of the coefficients. 

It should also be pointed out that even imperfect po- 
tentials could be quite useful for augmenting a fit of 
diffraction data. A combination of energy and diffrac- 
tion data could overcome minor spurious features in re- 
finements p^ ]. 

Our results strongly suggest that, in this alloy system, 
stability of the quasicrystal requires a ternary not only 
to tune the electron density, as in the Hume-Rothery pic- 
ture, but also because each species fills a particular type 
of site. Thus it would be highly desirable to explore some 
other composition regimes in simulations. The Co-rich 
"basic Co" subphase of d(AlNiCo) would be of particu- 
lar interest. It exhibits much stronger diffuse scattering 
than Basic Ni, particularly halfway between the layers of 
Bragg peaks in reciprocal space, indicating a local ten- 
dency to doubling the c periodicity (8 A structures). 

Our constraint that the atoms sit on ideal sites in the 
initial stage of energy optimization has considerable po- 
tential to distort the conclusions. Clearly, atom displace- 
ments have a very strong effect so it is crucial to include 
these for a final answer. Small adjustments of position 
can gain as much energy as swaps of species, and in some 
cases the optimal sites may lie halfway between ideal 
sites, or between ideal layers - the so-called "puckered" 
layers We do allow such relaxations in later stages, after 
determining a decoration, but it appears possible that 
a different decoration would be obtained, if we allowed 
relaxations in the early stages. This seems to be the 
most serious drawback of the lattice-gas approach 
we used. Its advantages are that the systematic Monte 
Carlo exploration is not biased by its users' prejudices, 
and furthermore it could readily be adapted to T > 
(the temperatures at which the quasicrystal phase is sta- 
ble). Nonzero temperature probably affects the result 
just as strongly as relaxed positions do. 

Within the framework of pair potentials and tile- 
decoration representation of the structure, a different ap- 
proach is possible [ jj"7| , whereby relaxed sites are used in 
the main discovery process. To implement this so that 
only valid tile decorations get explored, one performs re- 
laxations with the decoration-equivalent atoms of a given 
class ("orbit") being constrained to move together. How- 
ever, that method is ad-hoc, depending on its users' ed- 
ucated guesses as to the tilings and decorations to be 
tried. Future work should strive to blend the best fea- 
tures of these two approaches to using decorations with 
potentials for the discovery of structures. 
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FIG. 2. Two bilayers (thickness 2c), with periodic bound- 
ary conditions in all directions, after a tile flip in one bilayer, 
followed by molecular dynamics and relaxation to an energy 
minimum. Black lines are HBS supertile edges (length r 2 ao). 
The dashed circles are imperfect 20 A decagons - this cell is 
too small to contain complete ones. 
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FIG. 3. Probability density of Al atom centers during a 
molecular dynamics simulation of the same configuration as in 
Fig. [| TM locations are indicated by triangles, with the same 
coloring and size conventions as the previous figures. Edges of 
the large HBS tiles are marked by lines, (a) Projection onto 
xy plane, (b) Projection onto the yz plane of the narrow slab 
outlined by a dashed rectangle in (a). The pattern repeats 
with the z period 8 A. (Note this slice actually cuts through 
an H tile but not a 20 A cluster.) 



FIG. 1. Lowest-energy configuration obtained from a sim- 
ulation with variable occupancy interiors. Circles are col- 
ored with the following convention: empty=Al, gray=Co, 
black=Ni; the size of the symbol represents the z coordinate. 
Shaded small hexagon tiles, occupied by either NiNi pairs or 
AICo pairs, form the edges of an HBS supertiling. Dashed 
circle indicates an 20 A decagon cluster. This is the same 
configuration as Fig. 4(b) of Ref. Q . 
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